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We consider the interaction of electromagnetic radiation of arbitrary polarization 
with multi-level atoms in a self-consistent manner, taking into account both spatial 
and temporal dependencies of local fields. This is done by numerically solving the 
corresponding system of coupled Maxwell-Liouville equations for various geometries. 
In particular, we scrutinize linear optical properties of nanoscale atomic clusters, 
demonstrating the significant role played by collective effects and dephasing. It is 
shown that subwavelength atomic clusters exhibit two resonant modes, one of which 
is localized slightly below the atomic transition frequency of an individual atom, 
while the other is positioned considerably above it. As an initial exploration of future 
applications of this approach, the optical response of core-shell nanostructures, with 
a core consisting of silver and shell composed of resonant atoms, is examined. 

PACS numbers: 42.50. Ct, 78.67.-n 



I. INTRODUCTION 

Nanoscale optical materials have long been attracting considerable attention due to 
many important applications ranging from optical nanodevices [1], plasmonic circuitry [2], 
nanoscale sources of coherent radiation [3], single atom/molecule manipulation bio- 
medical applications ^ [6], and many others. Among such exciting applications lies the 
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yet-to-be explored sub-field of nanoscale optical atomic and molecular physics that deals 
with ensembles of atoms or molecules interacting with dielectric nanoparticles (NP) and 
their assemblies. Such systems are characterized by a significant spatial dependence of 
evanescent electromagnetic (EM) fields on the dielectric environment, providing the means 
to control the behavior of molecular systems by the combination of large EM fields and large 
field gradients, exemplified by recent work [714T2] on the EM field associated with metal NP 
dimers and their dependence on particle sizes and interparticle distance. At their core, these 
phenomena rely mostly on the excitations of surface plasmon-polariton (SPP) resonances 
[13] in systems comprising metal NPs and their arrays, [13] as well as other nanoscale metal 
surfaces, such as subwavelength diffraction gratings, whose optical properties depend sen- 
sitively on their surface topology and material parameters [15]. Their studies have lead to 
many applications such as coherent EM energy transport in space [16], surface enhanced 
Raman spectroscopy (SERS) [17] and tip-enhanced microscopy [T8]. Recent attention has 
focused on optical control scenarios, ranging from coupled exciton-plasmon dynamics in 
semiconductor nanodots [TUH27] and in molecular aggregates [2HHM] where metal NPs af- 
fect excitation energy transfer between molecules, to optical trapping of single atoms or 
molecules [371 - HO] . Such applications are facilitated by the possibility to control the geom- 
etry of nano-materials (NPs sizes, their relative arrangement, etc.) with an outstanding 
precision [IT] . 

While theoretical and computational methodologies for studying these phenomena have 
advanced considerably, the consequences of mutual feedback between molecular excitations 
and metallic SPP resonances are not well understood, especially when one probes systems 
comprising both metallic nanostructures and semiconductor or molecular particles or lay- 
ers. An often used simple description is based on assigning a dielectric response function 
to the semiconductor /molecular component and solve the electromagnetic problem for the 
corresponding composite dielectric. While such an approach can be useful for describing 
the effect of the molecular environment on the metallic plasmon, it cannot be used to de- 
scribe energy-transfer, relaxation, and spontaneous emission in the excitonic (molecular or 
semiconductor) system. It is hence important to develop a self-consistent description of 
the electromagnetic response of such systems. Such approach has to take into account the 
electrodynamics of the radiation field and the quantum dynamics of the molecular system 
in a self-consistent manner. This can be accomplished by solving simultaneously Maxwell's 
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equations for the radiation field and the Liouville equation for the molecular density matrix, 
including the molecular polarization current in the former and the molecules-field interaction 
in the latter. 

First attempts to consider numerically coupled Maxwell-Bloch equations have been ini- 
tiated by Ziolkowski et al. [l2l US] for simple two-level atoms in one and two dimensions 
utilizing finite-difference time-domain (FDTD) technique. Later on this approach has been 
extended to three dimensions [H] . Although these works contain interesting and important 
physics, they are limited to ensembles of two-level systems. Consideration of multilevel 
systems is critical for modeling of nano-lasing, which has to include at least three lev- 
els. Moreover, the proposed numerical implementation results in noticeably, long execution 
times. The scheme we propose is more efficient as discussed below. Similarly Neuhauser et 
al. proposed another approach where the authors coupled Maxwell's equations to the 
Schrodinger equation describing a molecule located in the closed proximity of a metal NP. 
These works, however, cannot in their present form include relaxation and dephasing effects, 
which, as we demonstrate, are very important. 

In this paper we describe a numerical implementation of such a model, using the method- 
ology developed by Ziolkowski et al. [l2l H3] as our starting point. This model captures 
collective effects that play pivotal role in electrodynamics of nano-systems, as well as the 
counterbalancing effect of dephasing processes. The questions to be addressed are: 

1. How does a size of the system affect scattering/absorption of EM radiation? 

2. What is a role of dephasing and relaxation effects? 

3. When does one observe collective response of atoms to external EM excitation? 

These and other closely related questions are not only important from the fundamental 
point of view (how optically induced interatomic or intermolecular interactions depend on 
structural/material parameters), but also essential for general understanding of optics of 
many-body systems. 

In this regard it should be pointed out that although the technique we propose in this 
paper is utilized to capture collective effects of quantum particles in the linear response 
regime, it can easily be applied to nonlinear systems. 

The paper is organized as follows: Section [IT] discusses our computational approach, based 



on coupled Maxwell-Liouville equations in the mean-filed approximation. In Section III we 
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provide details of the numerical implementation; Section |IV] describes and discusses the 
results of our numerical studies. Our main conclusions along with future research outlook 
are presented in Section |Vj 

II. MODEL 

We consider a general problem of a system of quantum particles (further referred to as 
atoms) interacting with EM radiation. We start from the time domain Maxwell's equations 
for the dynamics of the EM fields, E and H 

Ho-^ = -V X E, (la) 

^o^=Vx#-J, (lb) 

where /Xq and Eq are magnetic permeability and dielectric permittivity of free space, respec- 
tively. In spatial regions occupied by a metal nanostructure (such as metal NP, for instance) 



the equation ( lb ) is evaluated in the standard way from the metal dielectric dispersion |16] . 
In the present study the dispersion of dielectric constant of metal, b{u), is taken in the form 
of the Drude model 

= - (2) 

with numerical parameters describing silver for the wavelengths of interest Er = 8.26, Up = 
1.76 X 10^^ rad/sec, 7 = 3.08 x 10^^ rad/sec. The time evolution of the current density J in 
metal regions is then 

8 J 

— =aJ + bE, (3) 

where a = —7 and b = EoU^. 

In the spatial regions occupied by atoms, the mutual interaction between the atomic 
system and the EM field is accounted for in a self-consistent manner as follows: first, the 



current density in Eq. (lb) is expressed in terms of the macroscopic polarization of the 
atomic system, P (r , t) 

- dP 

The latter is given by 

P = na{fl), (5) 
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where 

(/i) = Tr(p/i) (6) 

is the expectation value of the atomic dipole moment and Ua is the atomic density. Eqs. 
Q - (|6]) constitutes the main approximation of the present approach, whereupon the local 
polarizability is expressed in terms of the local atomic density multiplied by the local aver- 
aged single atomic dipole. The time evolution of the latter is obtained from the evolution 
of the single atom density matrix (described below) in the presence of the EM field, thus 
providing a self-consistent description of the field-matter dynamics. 

Next consider the atomic subsystem. While our ultimate goal is to study realistic 3- 
dimensional systems, the present study focuses on nanoscale atomic clusters in two dimen- 
sions, taken to lie in the XY plane [IE]. The incident radiation field is represented by 
a transverse-electric (TE) mode with respect to ^-axis. It is characterized by two in-plane 
electric field components, and Ey, and one out-of-plane magnetic field component, H^. To 
account for the (two-dimensional) spherical symmetry of the atomic polarization response, 
the atoms are described as 3-level systems: an s-type ground state and two degenerate 
p-type excited states of px and Py character (as depicted in Fig. [l]A.). In anticipation of pos- 
sible generalizations to more complex models involving multilevel systems we use, in what 
follows, a basis of angular momentum wavefunctions with quantization axis in the z direc- 
tion, |1) = |s), |2) = {\px) +i\Py)) /v^, |3) = {\px) — 'i'\Py)) / V^i with optical transitions 
corresponding to AJ = ±1 and AM = ±1 selection rules. The corresponding Hamiltonian 
is 



^ VL^it) -n+{t)^ 



H = Ho-fI- E{t) 



(7) 



9.+ it) hhJa 

^ — r2_ (t) I 

where huja is the atomic energy transition, = figp {Ex (t) ± iEy (t)) / -\/6, and fXsp is s — p 
the matrix element of the dipole moment operator, that, in principle, can be taken from 
experiment or calculated using standard quantum chemistry packages. Its Cartesian com- 
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Note that the dipole moment operator (|8]) differs from the one used in [13] by the factor of 

Vs. 

The mean field approximation Eq. ([s]) make it possible to describe the atomic system in 
terms of the single atom density matrix p, which satisfies the Liouville equation 

dp 



ih 



dt 



[H, p] — ihVp. 



(9) 



Eq. ([9]) describes the time evolution of an atom interacting with the radiation field and 
subject to (assumed Markovian) relaxation processes described by the f operator, which 
is taken in the Lindblad form [47j. The diagonal elements of this operator correspond to 
excited states lifetimes, while nondiagonal elements account for dephasing effects. 

Eqs. ([1]) - (joj) describe the time evolution of the atomic system and radiation field 
in a self-consistent way. Note that in the single-atom Hamiltonian H given by Eq. ^ 
interatomic interactions are absent. Still, the dynamics described by Eqs. ([T]) - (|9]) is not 
that of independent atoms: the self-consistent scheme accounts for all interactions between 
atoms that are associated with their mutual interaction with the radiation field, including 
excitonic (energy transfer) interactions, which are all important for elucidating the overall 
system response. Note that, assuming that the EM field does not significantly vary within a 
volume occupied by a single atom, the spatial dependence of the density matrix in Eq. ^ 
depends parametrically on position via the EM field variables. 

Eqs. (|9])-([8]) lead to the following equations for the atomic density matrix elements 

dpii 



dt 



iu+ (pi2 + pIs) - iuj^ (pi3 + P12) + 7i (P22 + P33) 



dp 
~dt 



12 



i(^aPl2 - (p22 - Pll) + «W+P23 - 72^12, 



(10a) 
(10b) 
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= iUaPn + (p33 - pii) - iuj-p2z - 72P13, (10c) 

dp22 



dt 

dp23 

dt 

dp33 



-- iuj-pl2 - iuj+pi2 - 71P22, (lOd) 
-iu+ (pi3 + PI2)- 272P23, (lOe) 
= iu^Pis - i(^+P*i3 - 7iP33, (lOf) 



dt 

where oj± = Q±/h, 72 = 7p + 71/2 with 7p denoting the pure dephasing rate due to envi- 
ronmentally induced random fluctuations in the atomic energy spacing. As noted above, in 



Eq. (10) we denoted the ground state as |1) and the excited states | J = 1,M = —1) and 
\J = 1, M = 1) as |2) and |3), respectively. 



Finally, using Eqs. rtSl) and (10) we obtain the macroscopic polarization current (time 



derivative of Eq. rt5|), which enters Ampere's law (lb) 



^ - n ^ flla) 
dt dt ' 



where, back in cartesian coordinates 



- {p22 - P33) - [(Wa + ^72) (Pl2 - P13) - (Wa " «72) {pu " Pls)] ' (^b) 



dt 3h '^'^ ^'"^ V6 

^ = (P22 - P33) + ^ [{oOa + ^l2) (Pl2 + P13) + (C^a " ^72) (p^^ + P13)] > (^c) 

We end this section with two comments. First, as already pointed out, Eqs. (10) a mean 
fleld description of a system of atoms interacting with the EM fleld. In this approximation 
a single atom interacts with other atoms through the electromagnetic fleld associated with 
their mean local density. Obviously, such an approach cannot account for speciflc atom-atom 
correlations, but it can describe collective effects in a system of atoms resulting from their 
interaction with the EM fleld. 

Second, we note that this procedure can be easily generalized to yield the analogous 
3-dimensional coupled Maxwell-Liouville equations. To maintain spherical symmetry we 
would need to include an additional atomic level |J = 1,M = 0), which is coupled to the 
ground atomic state by Ez [44j. Obviously, it is also possible to expand the atomic basis 
and consider additional excited manifold starting with \ J = 2,M). Although the number 



of equations similar to Eqs. (10) grows signiflcantly, modern analytical computer packages 
such as Mathematica [19] can easily handle necessary algebra and subsequent computer 
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coding. For a molecule without rotational symmetry, average over angular distribution in 
the calculation of P from Eq. (j5| may be used when relevant. 

III. NUMERICAL APPROACH 

To solve the system ([T]), ( [lO| , (11) of coupled Maxwell-Liouville equations we employ a 



generalized FDTD technique (SD]. Within the FDTD, both the electric, E, and magnetic, H, 
fields are propagated in time and space by directly discretizing Maxwell's equations ([T]). This 
approach has several attractive technical features, including its numerical stability and the 
explicit description of the magnetic field. The latter is especially important if one considers 
structures with sharp corners, at which the tangential components of H have singularities. 
For the our purposes, the main advantage of the FDTD approach is that the boundary 
conditions (i.e. the continuity of the tangential H and E components) are automatically 
maintained at all grid points owing to the use of the Yee cell [51] . This allows straightforward 
programming of the complex geometries. 

For simulations of open systems, one needs to impose artificial absorbing boundaries in 
order to avoid numerical reflection of outgoing EM waves back to the simulation domain. 
Among the various approaches that address this numerical issue, the perfectly matched 
layers (PML) technique is considered to be the most powerful [52] ■ It reduces the reflection 
coefficient of outgoing waves at the simulation region boundary to ~ 10~^. In essence, the 
PML approach surrounds the simulation domain by thin layers of non-physical material that 
efficiently absorbs outgoing waves incident at any angle. We have implemented the most 
efficient and least memory variant of the method, the convolution perfectly matched layers 
(CPML) absorbing boundaries pBj. Through extensive numerical experimentation, we have 
empirically determined optimal parameters for the CPML boundaries that lead to almost 
no reflection of the outgoing EM waves at all incident angles. 

In the calculations reported below we consider structures with a characteristic size 
much smaller than the incident wavelength. Hence it is a good approximation to excite 
such systems using a plane wave. The latter is accomplished via implementation of total 
field/scattered fleld approach [5U] within the FDTD. 

We partition the FDTD scheme onto an array of parallel grid slices by dividing the cubic 
simulation cell into M xy slices, where M is the number of available processors. Point- 
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to-point message passing interface (MPI) communication subroutines [51] are implemented 
at the boundaries between slices. The number of xy planes in each slice usually varies in 
the range from 15 to 20. All simulations are performed on the home-build 128-core AMD 
Opteron based cluster at Arizona State University [55] . 

The numerical implementation of the proposed scheme is as follows 

1. In the spatial regions occupied by atoms the Maxwell equations are solved utilizing the 
standard FDTD algorithm. First, magnetic field is updated according to Faraday's 



law, Eq. (la). Next, using Ampere's law, Eq. (lb), we update the electric field with 
the macroscopic polarization current density, Eqs. (11), which is calculated using the 
density matrix of the previous time step. The EM fields in the regions occupied by 
metal are updated according to the auxiliary differential equation method, [HU] Eq. 
(§. 

2. With the knowledge of local electric field components (stored in memory at two pre- 
vious time steps) we update the density matrix at each spatial point on the grid 



according to Eq. (10) using the fourth order Runge-Kutta scheme 



3. Finally, with the knowledge of the electric field components and the updated den- 
sity matrix we calculate the macroscopic polarization current, at each grid point 



according to (11 ). 



We have verified this scheme using several test cases. A most important test of numerical 
stability is to check that the condition Tr(p) = 1 is maintained at each time step. In all 
simulations this condition was perfectly satisfied with almost no dependence on incident field 
amplitude and other physical parameters. Another interesting test was to demonstrate the 
absence of self-interactions in our calculation. Such interactions often appear spuriously in 
mean field calculations, whereupon a particle interacts with its own contribution to the mean 
density. In the present situation, however, the field produced by the oscillating dipole of a 
given atom propagates away from this atom and can affect it only through the polarization 
induced in other atoms or (in different setups) through reflection from the boundaries, both 
physically valid phenomena. We have verified that direct self-interaction is indeed absent 
in our calculation by solving Eqs. ([T]) and (10) for the case where the system occupies 



a single grid point: the same solution is obtained whether or not the polarization source 
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term is included in Eq. (lb). Finally, we have compared our results and execution time of 
the proposed integration scheme with those obtained by Ziolkowski et al. |43j. We have 
implemented the numerical approach based on the predict-corrector method and atomic 
basis as used in [13] and compared it with ours (keeping in mind that /i^p has to be re- 
normalized by the factor of -\/3). The simulation data obtained using both approaches were 
in excellent agreement. However execution times for the codes employing our approach were 
noticeably smaller. 

With the solution of Eqs. ([T]), ( p!o| , (11) obtained in this way, the following observables 
can be calculated: 

(1) The scattered radiation can be computed as the difference between the total and the 

incident EM fields. At any detection point, e.g. that depicted as a red diamond in Fig. 
[Tj3, we can calculate the Poynting vector components associated with the scattered 
EM field. This may be integrated over a spherical boundary surrounding the atomic 
system to yield the total scattered radiation. These calculations can be accomplished 
in a transient mode to give the time dependent response to an incident EM pulse, 
or in a steady state mode that yields the long time steady state response to CW 
incident radiation of a given frequency, cjj^. In the latter case we need to propagate 
the Maxwell-Liouville equations under the incident CW radiation until steady-state is 
reached. 

(2) Generally, for a given incident frequency the outgoing field may exhibit different 

frequencies c^out (with amplitudes obtained by Fourier-transforming the scattered sig- 
nal), making it possible to obtain the outgoing steady state flux (Poynting vector) in 
a given direction at any Wq^^ for a given incident frequency. Integrating the outgoing 
flux over Wout displaying the result as a function of yields the absorption 
spectrum of the atomic cluster. 

(3) A much easier way to calculate the steady state absorption lineshape is by calculating 

the steady state relaxation flux of the excited state populations according to 



dE 
~dt 



hwali j d^r (p22 (r) + P33 (r)) , (12) 



where the integral is taken over cluster volume [S^. Eq. (12) expresses for any given 



incoming frequency the rate of energy dissipation by the molecular system, which at 
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steady state should be equal to the rate of energy absorption at that frequency. 

(4) For some applications the short pulse method (SPM) [46j can save a substantial comput- 
ing effort. In this approach we use an ultra-short incident pulse with a wide bandwidth, 
which is almost flat in the spectral region of interest (in our simulations the incident 
pulse duration was set at 0.36 fs, which corresponds to flat spectrum throughout the 
frequency domain considered in this paper). Such pulse can be represented as a co- 
herent linear superposition of CW plane waves with different Wj^. We then propagate 
Maxwell-Liouville equations for several picoseconds (the total propagation time has to 
be significantly longer than the lifetime of excited states of an atom, I/71), and take 
the Fourier transform of the calculated field. Under conditions of linear response and 
elastic scattering (namely, when ^out ~ '-^in)' Fourier component at frequency a; 
contains all the information relevant to a CW process at frequency uj. 

When applicable, the SPM can save substantially on computation time, since it yields the 
system response at many frequencies from one short time computation. It is important to 
understand its shortcomings. Two limitations, mentioned above, are obvious: this method is 
applicable only in linear response and only when the light scattering process is elastic. Both 
limitations are associated with the requirement that a given incident frequency can give 
rise to response only at the same frequency. A third limitation is important in the present 



context, because Eqs. (10b) and (10c) are inherently nonlinear. Linearity is obtained when 
P22 and P33 may be neglected, and pn taken constant in these equations. In this case, the 
steady state solutions for pi2 and pi3, and therefore the polarization (Eq. ([s])), oscillate with 
the incident frequency as required. Care has to be exercised if one attempts to calculate 
the steady state excitation, P22 + P33, in this method. An approximation may be obtained 
if the incident frequency is close to atomic transition frequency, Ua, by taking the Fourier 
transform at w = 0, i.e. the time average, of the resulting time dependent signal. 

The method, as described, cannot describe spontaneous emission, i.e. fluorescence, since 
the latter is a quantum effect associated with the quantum nature of the radiation field. It 
has been demonstrated [58j that the effect of spontaneous emission can be accounted for 
partially by imposing a classical stochastic field on the system. Obviously, a classical EM 
noise cannot mimic vacuum fluctuations; in particular it can induce excitation of ground 
state molecules while vacuum fluctuations can lead only to radiative damping of excited 
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molecules. One can use this trick to study time evolution in a system which is initially 
inverted, that is, all molecules in the excited state. In this case, the induction of emission by 
the EM noise will soon lead to a dominant signal of induced emission that does not depend 
much on the nature of that noise, provided the latter is weak enough. This method has been 
succesfully used to simulate superradiance [58l 159] and gain |60] within the FDTD approach. 
However, this method cannot be used to generate fluorescence in a system of mostly ground 
state molecules (p22 + P33 ^ Pii), where the main effect of such noise will be to induce 
unphysical molecular excitation. 



IV. RESULTS AND DISCUSSION 

The simulation setup is shown in Fig. [T|3: an atomic cluster is excited with an x-polarized 
low intense plane wave propagating in the y direction. We use low intensity incoming fields 
(our incident electric field was fixed at 1 V/m) in order to insure linearity of the system 
response |6T]. Eqs. ([T]), ([To|, (11) are then evolved to yield the electric and magnetic fields 



as functions of time and position. In the current studies we focus on the y -component of 
the Poynting vector, 5*^ ~ E^Hz, as the observable of interest. 

Fig. [2] presents direct comparison of two methods: SPM approach and CW scheme. 
The scattering signals obtained from the two methods are in excellent agreement. Also 



shown is the absorption spectrum from the CW calculation using Eq. (12). The absorption 
lineshape (normalized so that it matches the scattering signals at the peak) also leads to 
nearly identical lineshape, except that it exhibits a more pronounced resonance near the 
atomic transition frequency. 

It is not surprising, but still providing a consistency check, that at the very low densities 
{ua < 10^^ m~^) and zero dephasing our simulations are in the perfect agreement with 
the Clausius-Mossotti approximation described above. Moreover, calculations, in which the 



polarization current term in (lb) was neglected (i.e. the atoms are not coupled to each 
other through their mutual interaction with the radiation field but rather driven only by 
external incident radiation), were in a perfect agreement with the data produced by the full 
self-consistent computations in the limit of low densities. 

Fig. [3] summarizes main results of our SPM calculations (note that we performed direct 
comparison of SPM data with that obtained via CW scheme for every set of parameters 
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discussed below). First the dependence of the scattering intensity on the density of atoms 
in the cluster is depicted in Fig. The atomic transition frequency is fixed at Ua = 3.1 eV 
(see figure caption for the rest of simulation parameters). The scattering radiation clearly 
exhibits two resonances; one, a relatively weak response, close to (slightly below) the atomic 
transition frequency. The other is a strong and broad peak at a higher frequency that moves 
to the blue at larger atomic densities. Additional simulations presented in Fig. |4] show that 
the intensity of the high frequency mode (unlike the low frequency one) scales as n"^ for the 
low density of atoms suggesting a possible collective nature of the peak. It should be noted 
that this collective mode is noticeably wider than the atomic transition resonance. 

The dependence of the scattering intensity on cluster's size is shown in Fig. [3^. The 
low energy resonance exhibits a red shift, when a radius of the cluster increases, while the 
resonant frequency of the high energy mode does not change with cluster's size. It, however, 
becomes significantly wider for larger clusters, which has been observed experimentally [62] . 

One of the advantages of the present calculations over the standard approach based on 
a dielectric model is the ability to examine the influence of the dephasing rate on optical 
properties. Fig. [3p shows simulation results obtained at three pure dephasing rates, jp, 
including the case without pure dephasing. We should note that at small 7p numerical 
simulations tend to become hard to converge at the frequencies near the collective resonance 
(in our case this occurs for 7^ < 4 x 10^^ s~^)- While for relatively high dephasing rates 
regular simulations require spatial steps, 6x, on the order of 1 nm, the case without pure 
dephasing should be explored at 6x < 0.2 nm. It is interesting to note that the collective 
mode, while decreasing its width with the decrease of the dephasing rate, is still significantly 
wider than the atomic transition peak. In contrast, the low frequency peak becomes narrower 
with decreasing dephasing rate, with its width approaching 71. Note that the scattering 
actually shows a dip at the atomic frequency, with the low frequency peak slightly below it. 

Fig. [3p explores how the scattering is affected by the matrix element of the atomic 
dipole moment, fisp- Not surprisingly, the result for increasing fisp is quahtatively similar 
to that obtained with increasing atomic density. It is seen that larger /i^p results in blue 
shift of the higher frequency collective mode, and in a red shift of the lower frequency mode. 
The former shows a quadratic dependence of the resonant frequency on the dipole moment, 
which has been theoretically discussed in the case of a sphere with uniformly distributed 
linear quantum dipoles [55] . 
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It is useful, for the sake of comparison, to consider the simplest theoretical description for 
the optical response of our system, by modeling it as a dielectric particle with a dielectric 
response function taken from the Clausius-Mossotti expression 

l + 2x 



I — X 

An „ 



(13) 



where x = ^na, n is the number density of atoms, and a is the atomic polarizability. The 
absorption hneshape can be calculated [53] as the ratio between the dissipated power, Pdiss = 
(1/2) J d^ra\Ef and the incident flux, Jj^ = c|i?j^|^/ (Svr), where a = (wj^/ (47r)) lm{e) 
is the conductivity, e is the dielectric response function, E is the electric field in the particle, 
E^j^ is the incident electric field and c is the speed of light. For a small spherical particle 

m 

E = ^—E;^. (14) 

Using these expressions we obtain the absorption cross-section of a small spherical particle 
of volume Vt in the form 



_ J^diss _ Q ^m 
a — J — ^ ' 

"'in 



Im (e) = ^^^nVt Im (a) . (15) 



e + 2 

The imaginary part of the molecular optical polarizability is essentially a Lorentzian reso- 
nance peaked at the atomic transition frequency uOa- We see that the absorption cross-section 
in the Clausius-Mossotti approximation is proportional to the number of particles and to 
the absorption of a single particle, as would be predicted for a system of non-interacting 
particles. On the same level of theory, the dipole induced on the particle is [65] 



/i = — Tl^^^P^'wy — ^^«-^in- (16) 



2 

And, since the scattered light is proportional to |/i| , it is predicted to go like the square of the 
number of particles, and to have a similar resonant behavior as a function of the incident 



frequency. Eqs. (15) and (16) describe essentially a system of non-interacting particles, 
occupying a volume with linear dimensions much smaller than the radiation wavelength, 
that respond coherently to the incident radiation. This approximation is valid at low atomic 
density. We will see below how dephasing and through-field interatomic interactions at 
higher densities affect this behavior. 

While the calculation procedure applied here provides a route to explore the effect of 
dephasing on the optical response of atomic and molecular clusters, further studies will be 
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needed in order to determine how much of this effect was indeed captured in our simulations. 
The destruction of phase is affected within our calculations on the mean field evolution of a 
single atom, feeling the effects of others through their mutual interaction with the radiation 
field. It is not obvious that this mean field implementation can capture the full physics of 
atoms going out of phase from each other. We leave this important technical question to a 
future study. 

It is informative to explore the spatial dependence of EM intensity, / ~ -E"^ + E"^, at reso- 
nant conditions. Fig. |5] shows intensity distributions calculated using steady-state solutions 
of Maxwell-Liouville equations for the two resonance modes. Clearly the EM intensity at the 
lower resonant frequency is mainly localized on the surface of the cluster exhibiting dipole 
radiation pattern similar to EM intensity distributions seen at the plasmon resonance for 
a single metal nanoparticle fL3\. The collective high frequency mode is distributed in the 
entire volume of the particle, where all atoms coherently participate in the radiation pro- 
cess. This suggests that the high frequency mode is more collective in nature than the low 
frequency one, consistent with its large density dependent shift from the atomic frequency 
and its scaling at low densities. The low frequency mode, involving fewer surface atoms, 
may be more atomic in nature. It is important to note however that this is not a single 
atom response since it clearly shifts with increasing atomic density. 

Optical properties of molecular aggregates resonantly coupled to plasmonic materials 
have been a subject of extensive research for the past several years. Sugawara et al. [UU] 
demonstrated significant modification of transmission and reflection spectra of a gold film 
with deposited J-aggregates. It has been shown experimentally that SPP resonances no- 
tably affect molecular electronic structure leading to resonance splitting [33j. The latter 
was proposed to be used for controlling optics of such hybrid material using femtosecond 
laser pulses Moreover core-shell metal NPs with a shell comprised of optically active 
molecules have been recently studied experimentally [SB]- To demonstrate the generality of 
our approach we present simulations of the core-shell particle schematically depicted in the 
inset of Fig. |6j Here a silver nanoparticle is shelled by a resonant atomic layer, with atomic 
transition frequency equaled to the SPP resonance of silver, Ua = 3.61 eV. The optical prop- 
erties of silver are described within the Drude model with parameters as in [18]. The hollow 
atomic shell exhibits doubled collective mode (red dashed line in Fig. [6]), which corresponds 
to the symmetry of the problem and can be understood within the plasmon hybridization 
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model proposed for noble metal core-shell particles [69]. The important observation is a 
clear splitting of the SPP mode with additional strong peak centered near atomic transition 
frequency (blue dash-dotted line in Fig. |6|. The observed splitting as indicated in is 
due to the strong optical coupling of atoms with the SPP mode. 



V. CONCLUSION 



We have presented a self-consistent electrodynamical model based on coupled Maxwell- 
Liouville equations that takes into account arbitrary polarization of the incident field. The 
proposed model is applied to investigate linear optical response of nanoscale atomic clusters 
in two dimensions merging classical electrodynamics with quantum mechanical description of 
atoms. The calculations can capture collective effects that play pivotal role in electrodynam- 
ics of nano-systems and, with limitations discussed below, includes the effect of dephasing 
on the optical response of these systems. 

We have found that spherical atomic clusters exhibit two well-distinguished resonances. 
The low energy resonance is close to the atomic transition frequency of individual atoms. 
The EM intensity distribution at this resonance is localized near the surface of a cluster. The 
high energy mode, where all atoms in the cluster coherently participate in the scattering, 
has clear collective nature. The dependence of the scattering intensity on various parameters 
was considered. It was demonstrated that the pure dephasing plays an important role in the 
scattering and absorption. Moreover we successfully applied our formalism to more complex 
systems, which comprise a resonant atomic shell and a silver core. 

Applications of the proposed scheme are many and vast. They range from complete 
three-dimensional description of nanoparticles resonantly coupled to ensembles of quantum 
particles, to nonlinear optical phenomena at the nanoscale. We note that our scheme can 
be extended to molecular systems, where one may investigate Raman processes. At the 
same time, we have indicated physically significant open technical issues. One is the limited 
ability of an approach based on classical electrodynamics to describe spontaneous emission 
and hence fluorescence. To account for such phenomena we need to modify the Maxwell- 
Liouville equations, Eqs. ([T]), ( p!o| , (11), so as to take into account the quantum nature 
of the radiation field, possibly using the quantization schemes described in Refs. [70H72] 
or [731 El] (which were shown to be equivalent [75]). Another is the need to examine the 
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adequacy of mean field calculations of dephasing. All these will be subjects of our continuing 
studies. 
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FIG. 1: (Color online). Panel A shows the energy level diagram of a two-level two-dimensional 
atom with black arrows indicating optically induced transitions by the TE mode and red arrows 
representing spontaneous decay. Panel B depicts schematics of the simulations with the detection 
point shown as a red diamond (in the lower left corner) , where the y-component of the Poynting 
vector is calculated. 
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FIG. 2: Scattering intensity as a function of the incoming frequency uj^^ calculated within SPM 
approach (solid line) and CW scheme (circles). Normalized absorption (see Eq.(|12[)) is shown as 
squares. Simulations are performed for the cluster with the following parameters: Ua = 3.1 eV, 
= 25 nm, = 7 X 10^5 m'^, 71 = lO^^ s'^, jp = 10^^ g-i^ = 25 Debye. 
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FIG. 3: (Color online). Linear optics of atomic clusters with uja = 3.1 eV, 71 = 10^^ s~^. Panel 
A: scattering intensity as a function of the incident frequency for three atomic densities, - black 
solid line: = 2.5 x 10^^ m~'^ (ideal gas at atmospheric pressure and room temperature), red 
dashed line: = 7 x 10^^ m~^, and blue dash-dotted line: Ua = 10^^ m^'^. Other parameters are: 
R = 25 nm, = 10^'^ s"^, fisp = 25 Debye. Panel B: same is in panel A for three cluster's radii, 
R - black (solid), red (dashed), and blue (dash-dotted) lines show scattering intensity results for 
i? = 15 nm, 25 nm, and 35 nm, respectively. Other parameters are: = m ^, 7p = 10^^ s ^, 
Hsp = 25 Debye. Panel C: same as in panels A and B (now shown in logarithmic scale), but for 
four pure dephasing rates, 7^ - black (solid) line: 7^ = 2 x 10^^ s~^, red (dashed) line: l7p = 10^^ 
s~^, blue (dash-dotted) line 7^ = 2 x 10^^ g-i^ and green (dotted) line shows the data without pure 
dephasing 7p = s~^. Other parameters are: R = 25 nm, = 10^^ m~^, fisp = 25 Debye. Panel 
D: same as in panels A-C, but for three values of the matrix element of the dipole moment, figp ~ 
black (solid) , red (dashed) , and blue (dash-dotted) lines correspond to atomic systems characterized 
by fJ-sp = 10 Debye, fisp = 25 Debye, and /i^p = 40 Debye, respectively. Other parameters are: 
R = 25 nm, = 10^^ m-^, = 10^2 g-^, 7p = 10^^ s"^ 
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FIG. 4: (Color online). Scattering intensity (shown in double logarithmic scale) at the high fre- 
quency resonance, u)res, as a function of the atomic density, n^, for three pure dephasing rates, jp 
- blue circles: 7^ = 2 x 10^^ s~^, black rectangles: jp = 10^^ s~^, and red triangles: 7^ = 2 x 10^^ 
s~^. Other parameters are: coa = 3.1 eV, 71 = 10^^ s^^, R = 25 nm, tIq = 10^^ m~^, fi^p = 25 
Debye. The dashed straight lines represent fitting for each set of data demonstrating nearly ideal 
quadratic dependence on at low densities. 
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CJOin = 3.09 eV CJOin = 3.18 eV 




FIG. 5: (Color online). Spatial distributions of EM intensity (normalized to the incident intensity) 
in logarithmic scale at the low frequency resonance (left panel) and the high frequency resonance 
(right panel). Cluster's parameters are: Ua = 3.1 eV, 71 = 10"*^^ s^"*^, 7^ = 10^^ s^^, R = 25 nm, 
na = 10^^ m~'^, fisp = 25 Debye. 
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FIG. 6: (Color online). Scattering intensity as a function of the incident frequency for the core- 
shell particle shown in the inset. Black solid line shows the data for the silver nanoparticle without 
atomic shell, red dashed line presents simulations for the hollow atomic shell, and blue dash-dotted 
line demonstrates results obtained for atomic shell with a silver core. Parameters of the simulations 



are: n„ = 2.5 x 10^^ m 
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